using System;
using Xunit;

namespace NPack.Tests
{
    ///<summary>
    /// Various tests involving computing known mathematical quantities 
    /// using the basic <see cref="DoubleDoubleComponent"/> arithmetic operations.
    ///</summary>
    public class DoubleDoubleComputeTest
    {

        [Fact]
        public void TestEByTaylorSeries()
        {
            DoubleDoubleComponent testE = ComputeEByTaylorSeries();
            double err = Math.Abs((double)testE.Subtract(DoubleDoubleComponent.E));
            Console.WriteLine(err);
            Assert.True(err < 64 * DoubleDoubleComponent.EPS);
        }

        ///<summary>
        /// Uses Taylor series to compute e
        /// e = 1 + 1 + 1/2! + 1/3! + 1/4! + ...
        ///</summary>
        /// <returns>an approximation to e</returns>
        private static DoubleDoubleComponent ComputeEByTaylorSeries()
        {
            DoubleDoubleComponent s = new DoubleDoubleComponent(2.0);
            DoubleDoubleComponent t = new DoubleDoubleComponent(1.0);
            double n = 1.0;
            int i = 0;

            while ((Double)t > DoubleDoubleComponent.EPS)
            {
                i++;
                n += 1.0;
                t = t.Divide(new DoubleDoubleComponent(n));
                s = s.Add(t);
                Console.WriteLine(s);
            }
            return s;
        }

        [Fact]
        public void TestPiByMachin()
        {
            DoubleDoubleComponent testE = computePiByMachin();
            double err = Math.Abs((double)testE.Subtract(DoubleDoubleComponent.PI));
            Console.WriteLine(err);
            Assert.True(err < 8 * DoubleDoubleComponent.EPS);
        }


        ///<summary>
        /// Uses Machin's arctangent formula to compute Pi:
        /// Pi / 4  =  4 arctan(1/5) - arctan(1/239)
        ///</summary>
        /// <returns>an approximation to Pi</returns>
        private static DoubleDoubleComponent computePiByMachin()
        {
            DoubleDoubleComponent t1 = new DoubleDoubleComponent(1.0).Divide(new DoubleDoubleComponent(5.0));
            DoubleDoubleComponent t2 = new DoubleDoubleComponent(1.0).Divide(new DoubleDoubleComponent(239.0));

            DoubleDoubleComponent pi4 = (new DoubleDoubleComponent(4.0)
                .Multiply(arctan(t1)))
                .Subtract(arctan(t2));
            DoubleDoubleComponent pi = new DoubleDoubleComponent(4.0).Multiply(pi4);
            Console.WriteLine(pi);
            return pi;
        }

        ///<summary>
        /// Computes the arctangent based on the Taylor series expansion
        /// arctan(x) = x - x^3 / 3 + x^5 / 5 - x^7 / 7 + ...
        ///</summary>
        ///<param name="x">the argument</param>
        /// <returns>an approximation to the arctangent of the input</returns>
        private static DoubleDoubleComponent arctan(DoubleDoubleComponent x)
        {
            DoubleDoubleComponent t = x;
            DoubleDoubleComponent t2 = t.Multiply(t);
            DoubleDoubleComponent at = new DoubleDoubleComponent(0.0);
            DoubleDoubleComponent two = new DoubleDoubleComponent(2.0);
            int k = 0;
            DoubleDoubleComponent d = new DoubleDoubleComponent(1.0);
            int sign = 1;
            while ((double)t > DoubleDoubleComponent.EPS)
            {
                k++;
                if (sign < 0)
                    at = at.Subtract(t.Divide(d));
                else
                    at = at.Add(t.Divide(d));

                d = d.Add(two);
                t = t.Multiply(t2);
                sign = -sign;
            }
            Console.WriteLine("Computed: " + at
                    + "    Math.atan = " + Math.Atan((double)x));
            return at;
        }

    }

    ///<summary>
    /// Tests basic arithmetic operations for <see cref="DoubleDoubleComponent"/>s.
    ///</summary>
    public class DoubleDoubleBasicTest
    {
        [Fact]
        public void TestNaN()
        {
            Assert.True(new DoubleDoubleComponent(1).Divide(new DoubleDoubleComponent(0)).IsNaN());
            Assert.True(new DoubleDoubleComponent(1).Multiply(DoubleDoubleComponent.NaN).IsNaN());
        }


        [Fact]
        public void TestAddMult2()
        {
            CheckAddMult2(new DoubleDoubleComponent(3));
            CheckAddMult2(DoubleDoubleComponent.PI);
        }

        [Fact]
        public void TestMultiplyDivide()
        {
            CheckMultiplyDivide(DoubleDoubleComponent.PI, DoubleDoubleComponent.E, 1e-30);
            CheckMultiplyDivide(DoubleDoubleComponent.TWO_PI, DoubleDoubleComponent.E, 1e-30);
            CheckMultiplyDivide(DoubleDoubleComponent.PI_2, DoubleDoubleComponent.E, 1e-30);
            CheckMultiplyDivide(new DoubleDoubleComponent(39.4), new DoubleDoubleComponent(10), 1e-30);
        }

        [Fact]
        public void TestDivideMultiply()
        {
            CheckDivideMultiply(DoubleDoubleComponent.PI, DoubleDoubleComponent.E, 1e-30);
            CheckDivideMultiply(new DoubleDoubleComponent(39.4), new DoubleDoubleComponent(10), 1e-30);
        }

        [Fact]
        public void TestSqrt()
        {
            // the appropriate error bound is determined empirically
            CheckSqrt(DoubleDoubleComponent.PI, 1e-30);
            CheckSqrt(DoubleDoubleComponent.E, 1e-30);
            CheckSqrt(new DoubleDoubleComponent(999.0), 1e-28);
        }

        private static void CheckSqrt(DoubleDoubleComponent x, double errBound)
        {
            DoubleDoubleComponent sqrt = x.Sqrt();
            DoubleDoubleComponent x2 = sqrt.Multiply(sqrt);
            CheckErrorBound("Sqrt", x, x2, errBound);
        }

        /*
      [Fact() ][Ignore]
      public void testTrunc()
      {
        checkTrunc(new DoubleDoubleComponent(1e16).Subtract(new DoubleDoubleComponent(1)),
                new DoubleDoubleComponent(1e16).Subtract(new DoubleDoubleComponent(1)));
        // the appropriate error bound is determined empirically
        checkTrunc(DoubleDoubleComponent.PI, new DoubleDoubleComponent(3));
        checkTrunc(new DoubleDoubleComponent(999.999), new DoubleDoubleComponent(999));
  	
        checkTrunc(DoubleDoubleComponent.E.Negative(), new DoubleDoubleComponent(-2));
        checkTrunc(new DoubleDoubleComponent(-999.999), new DoubleDoubleComponent(-999));
      }

      private void checkTrunc(DoubleDoubleComponent x, DoubleDoubleComponent expected)
      {
        DoubleDoubleComponent trunc = x.trunc();
        boolean isEqual = trunc.equals(expected);
        assertTrue(isEqual);
      }
        */
        [Fact]
        public void TestPow()
        {
            CheckPow(0, 3, 16 * DoubleDoubleComponent.EPS);
            CheckPow(14, 3, 16 * DoubleDoubleComponent.EPS);
            CheckPow(3, -5, 16 * DoubleDoubleComponent.EPS);
            CheckPow(-3, 5, 16 * DoubleDoubleComponent.EPS);
            CheckPow(-3, -5, 16 * DoubleDoubleComponent.EPS);
            CheckPow(0.12345, -5, 1e5 * DoubleDoubleComponent.EPS);
        }
        [Fact]
        public void TestReciprocal()
        {
            // error bounds are chosen to be "close enough" (i.e. heuristically)

            // for some reason many reciprocals are exact
            CheckReciprocal(3.0, 0);
            CheckReciprocal(99.0, 1e-29);
            CheckReciprocal(999.0, 0);
            CheckReciprocal(314159269.0, 0);
        }

        [Fact]
        public void TestBinom()
        {
            CheckBinomialSquare(100.0, 1.0);
            CheckBinomialSquare(1000.0, 1.0);
            CheckBinomialSquare(10000.0, 1.0);
            CheckBinomialSquare(100000.0, 1.0);
            CheckBinomialSquare(1000000.0, 1.0);
            CheckBinomialSquare(1e8, 1.0);
            CheckBinomialSquare(1e10, 1.0);
            CheckBinomialSquare(1e14, 1.0);
            // Following call will fail, because it requires 32 digits of precision
            //  	checkBinomialSquare(1e16, 1.0);

            CheckBinomialSquare(1e14, 291.0);
            CheckBinomialSquare(5e14, 291.0);
            CheckBinomialSquare(5e14, 345291.0);
        }

        private static void CheckAddMult2(DoubleDoubleComponent dd)
        {
            DoubleDoubleComponent sum = dd.Add(dd);
            DoubleDoubleComponent prod = dd.Multiply(new DoubleDoubleComponent(2.0));
            CheckErrorBound("AddMult2", sum, prod, 0.0);
        }

        private static void CheckMultiplyDivide(DoubleDoubleComponent a, DoubleDoubleComponent b, double errBound)
        {
            DoubleDoubleComponent a2 = a.Multiply(b).Divide(b);
            CheckErrorBound("MultiplyDivide", a, a2, errBound);
        }

        private static void CheckDivideMultiply(DoubleDoubleComponent a, DoubleDoubleComponent b, double errBound)
        {
            DoubleDoubleComponent a2 = a.Divide(b).Multiply(b);
            CheckErrorBound("DivideMultiply", a, a2, errBound);
        }

        private static void CheckErrorBound(String tag, DoubleDoubleComponent x, DoubleDoubleComponent y, double errBound)
        {
            DoubleDoubleComponent err = x.Subtract(y).Abs();
            Console.WriteLine(tag + " err=" + err);
            Boolean isWithinEps = (double)err <= errBound;
            Assert.True(isWithinEps);
        }

        /**
         * Computes (a+b)^2 in two different ways and compares the result.
         * For correct results, a and b should be integers.
         * 
         * @param a
         * @param b
         */

        static void CheckBinomialSquare(double a, double b)
        {
            // binomial square
            DoubleDoubleComponent add = new DoubleDoubleComponent(a);
            DoubleDoubleComponent bdd = new DoubleDoubleComponent(b);
            DoubleDoubleComponent aPlusb = add.Add(bdd);
            DoubleDoubleComponent abSq = aPlusb.Multiply(aPlusb);
            //  	System.out.println("(a+b)^2 = " + abSq);

            // expansion
            DoubleDoubleComponent a2dd = add.Multiply(add);
            DoubleDoubleComponent b2dd = bdd.Multiply(bdd);
            DoubleDoubleComponent ab = add.Multiply(bdd);
            DoubleDoubleComponent sum = b2dd.Add(ab).Add(ab);

            //  	System.out.println("2ab+b^2 = " + sum);

            DoubleDoubleComponent diff = abSq.Subtract(a2dd);
            //  	System.out.println("(a+b)^2 - a^2 = " + diff);

            DoubleDoubleComponent delta = diff.Subtract(sum);

            Console.WriteLine();
            Console.WriteLine("A = " + a + ", B = " + b);
            Console.WriteLine("[DoubleDoubleComponent] 2ab+b^2 = " + sum
                    + "   (a+b)^2 - a^2 = " + diff
                    + "   delta = " + delta);
            PrintBinomialSquareDouble(a, b);

            Boolean isSame = diff.Equals(sum);
            Assert.True(isSame);
            Boolean isDeltaZero = delta.IsZero();
            Assert.True(isDeltaZero);
        }

        static void PrintBinomialSquareDouble(double a, double b)
        {
            double sum = 2 * a * b + b * b;
            double diff = (a + b) * (a + b) - a * a;
            Console.WriteLine("[double] 2ab+b^2= " + sum
                    + "   (a+b)^2-a^2= " + diff
                    + "   delta= " + (sum - diff));
        }

        [Fact]
        public void TestBinomial2()
        {
            CheckBinomial2(100.0, 1.0);
            CheckBinomial2(1000.0, 1.0);
            CheckBinomial2(10000.0, 1.0);
            CheckBinomial2(100000.0, 1.0);
            CheckBinomial2(1000000.0, 1.0);
            CheckBinomial2(1e8, 1.0);
            CheckBinomial2(1e10, 1.0);
            CheckBinomial2(1e14, 1.0);

            CheckBinomial2(1e14, 291.0);

            CheckBinomial2(5e14, 291.0);
            CheckBinomial2(5e14, 345291.0);
        }

        static void CheckBinomial2(double a, double b)
        {
            // binomial product
            DoubleDoubleComponent add = new DoubleDoubleComponent(a);
            DoubleDoubleComponent bdd = new DoubleDoubleComponent(b);
            DoubleDoubleComponent aPlusb = add.Add(bdd);
            DoubleDoubleComponent aSubb = add.Subtract(bdd);
            DoubleDoubleComponent abProd = aPlusb.Multiply(aSubb);
            //  	System.out.println("(a+b)^2 = " + abSq);

            // expansion
            DoubleDoubleComponent a2dd = add.Multiply(add);
            DoubleDoubleComponent b2dd = bdd.Multiply(bdd);

            //  	System.out.println("2ab+b^2 = " + sum);

            // this should equal b^2
            DoubleDoubleComponent diff = abProd.Subtract(a2dd).Negative();
            //  	System.out.println("(a+b)^2 - a^2 = " + diff);

            DoubleDoubleComponent delta = diff.Subtract(b2dd);

            Console.WriteLine();
            Console.WriteLine("A = " + a + ", B = " + b);
            Console.WriteLine("[DoubleDoubleComponent] (a+b)(a-b) = " + abProd
                              + "   -((a^2 - b^2) - a^2) = " + diff
                              + "   delta = " + delta);
            //  	printBinomialSquareDouble(a,b);

            Boolean isSame = diff.Equals(b2dd);
            Assert.True(isSame);
            Boolean isDeltaZero = delta.IsZero();
            Assert.True(isDeltaZero);
        }


        private static void CheckReciprocal(double x, double errBound)
        {
            DoubleDoubleComponent xdd = new DoubleDoubleComponent(x);
            DoubleDoubleComponent rr = xdd.Reciprocal().Reciprocal();

            double err = (double)(xdd.Subtract(rr));

            Console.WriteLine("DD Recip = " + xdd
                              + " DD delta= " + err
                              + " double recip delta= " + (x - 1.0 / (1.0 / x)));

            Assert.True(err <= errBound);
        }

        private static void CheckPow(double x, int exp, double errBound)
        {
            DoubleDoubleComponent xdd = new DoubleDoubleComponent(x);
            DoubleDoubleComponent pow = xdd.Power(exp);
            Console.WriteLine("Pow(" + x + ", " + exp + ") = " + pow);
            DoubleDoubleComponent pow2 = SlowPow(xdd, exp);

            double err = (double)pow.Subtract(pow2);

            Boolean isOK = err < errBound;
            if (!isOK)
                Console.WriteLine("Test slowpow value " + pow2);

            Assert.True(err <= errBound);
        }

        private static DoubleDoubleComponent SlowPow(DoubleDoubleComponent x, int exp)
        {
            if (exp == 0)
                return new DoubleDoubleComponent(1.0);

            int n = Math.Abs(exp);
            // MD - could use binary exponentiation for better precision & speed
            DoubleDoubleComponent pow = new DoubleDoubleComponent(x);
            for (int i = 1; i < n; i++)
            {
                pow = pow.Multiply(x);
            }
            if (exp < 0)
            {
                return pow.Reciprocal();
            }
            return pow;
        }
    }
}